Mega-analysis of the Interoceptive Accuracy Scale (IAS)

Study 1

Code
library(tidyverse)
library(easystats)
library(patchwork)
library(lavaan)
library(ggraph)
library(tidySEM)
library(EGAnet)

Data Preparation

  • Murphy (2020)
    • https://osf.io/3m5nh
Code
df1a <- haven::read_sav("../data/Murphy2020/Study 1.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))

df1b <- haven::read_sav("../data/Murphy2020/Study 6 IAS.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))
  • Gaggero (2021)
    • https://osf.io/5x9sg
Code
load("../data/Gaggero2021/data.rda")
df2 <- data |> 
  mutate(Gender = as.character(gender)) |> 
  mutate(across(starts_with("IAS "), as.numeric)) |>
  rename(
    Age=age, Heart = `IAS 1`, Hungry = `IAS 2`, Breathing = `IAS 3`, Thirsty = `IAS 4`,
    Urinate = `IAS 5`, Defecate = `IAS 6`, Taste = `IAS 7`, Vomit = `IAS 8`, Sneeze = `IAS 9`,
    Cough = `IAS 10`, Temp = `IAS 11`, Sex_arousal = `IAS 12`, Wind = `IAS 13`, Burp = `IAS 14`,
    Muscles = `IAS 15`, Bruise = `IAS 16`, Pain = `IAS 17`, Blood_Sugar = `IAS 18`,
    Affective_touch = `IAS 19`, Tickle = `IAS 20`, Itch = `IAS 21`)
rm(data)
  • Campos (2022) - Study 1
    • https://osf.io/j6ef3
Code
df3 <- haven::read_sav("../data/Campos2022/Dataset_Test.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Sex == 1, "Male", ifelse(Sex == 0, "Female", "Other")))) |>
  rename(
    Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
    Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
    Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
    Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
    Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
  )
  • Todd (2022)
    • https://osf.io/ms354/
Code
df4 <- haven::read_sav("../data/Todd2022/CompleteData_new.sav") |>
  rename_with(\(x) str_remove(x, "IAS_"), .cols = starts_with("IAS_")) |>
  rename(
    Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
    Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
    Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
    Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
    Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
  )
  • Arslanova (2022)
    • https://osf.io/mp3cy/
    • note: Final.xlsx includes information about which participants passed the attention checks
Code
df5_attention <- readxl::read_excel("../data/Arslanova2022/Participants2.xlsx", 
    sheet = "FINAL") |>
  filter(ActorHR == 1) |>
  select(Subj.ID, Gender, Age) |>
  mutate(Gender = as.numeric(Gender), Age = as.numeric(Age)) |>
  filter(!is.na(Age)) |> #error on the documentation of this participant - attention check failed but not noted
  mutate(Gender = case_when(
    Gender== 0 ~ "Male",
    Gender== 1 ~ "Female" # based on paper reporting 65 women
  ))
New names:
• `Why No?` -> `Why No?...4`
• `Why No?` -> `Why No?...6`
• `` -> `...12`
• `` -> `...18`
• `` -> `...19`
• `` -> `...20`
• `` -> `...21`
• `` -> `...22`
• `` -> `...23`
• `` -> `...24`
• `` -> `...25`
• `` -> `...26`
• `` -> `...27`
• `` -> `...28`
Code
df5 <- readxl::read_excel("../data/Arslanova2022/Murphy_Iacc_scale.xlsx")  |>
  filter(str_detect(Question.Key, "IAC_")) |> 
  filter(str_detect(Question.Key, "-quantised")) |> 
  pivot_wider(names_from = Question.Key, values_from = Response) |>
  rename_with(\(x) str_remove(x, "IAC_"), .cols = starts_with("IAC_")) |> 
  rename_with(\(x) str_remove(x, "-quantised"), .cols = everything()) |> 
  rename(Sex_arousal = SexuallyAroused, Itch = Itchy, Sex_arousal = SexuallyAroused,
         Temp = HotCold, Tickle = Ticklish, Breathing= BreatheFast, 
         Affective_touch= PleasantAffectionate, Blood_Sugar= BloodSugar, 
         Muscles=TiredMuscles, Heart= FastHeart, Taste=Tastes) |> 
  rename(Subj.ID = "Participant.Public.ID") |>
  select(-starts_with("Participant")) 

df5 <- df5 |>
  filter(Subj.ID %in% df5_attention$Subj.ID) |>
  select(-Subj.ID) |>
  mutate(across(everything(), as.numeric))
  • Brand (2022)
    • https://osf.io/xwz6g/
Code
df6 <- haven::read_sav("../data/Brand2022/LatentVariableApproach.sav") |> 
  select(-SERIAL) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )
  • Brand (2023)
    • https://osf.io/3f2h6/
Code
df7a <- haven::read_sav("../data/Brand2023/Data_Giessen.sav") |> 
  select(-COUNTRY_OTHER) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(SEX == 1, "Male", ifelse(SEX == 2, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

df7b <- haven::read_sav("../data/Brand2023/Data_Mainz.sav") |> 
  select(-COUNTRY_OTHER, -EDUCATION_OTHER, -Sample) |> 
  mutate_all(as.numeric) |> 
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )
  
  
df7c <- haven::read_sav("../data/Brand2023/Data_PotVie.sav") |> 
  select(-ID) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IA02_01, Hungry = IA02_02, Breathing = IA02_03, Thirsty = IA02_04,
    Urinate = IA02_05, Defecate = IA02_06, Taste = IA02_07, Vomit = IA02_08, Sneeze = IA02_09,
    Cough = IA02_10, Temp = IA02_11, Sex_arousal = IA02_12, Wind = IA02_13, Burp = IA02_14,
    Muscles = IA02_15, Bruise = IA02_16, Pain = IA02_17, Blood_Sugar = IA02_18,
    Affective_touch = IA02_19, Tickle = IA02_20, Itch = IA02_21
  )
  • Lin (2023)
    • https://osf.io/3eztd
    • Note: No tickle because same Chinese character
Code
df8a <- haven::read_sav("../data/Lin2023/Study 1 & 3.sav") |>
  select(-sex) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(sex_dummy == 1, "Male", ifelse(sex_dummy == 0, "Female", "Other")))) |>
  rename(
    Age = age,
    Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
    Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
    Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
    Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
    Affective_touch = Touch, Itch = ITCH
  ) 


df8b <- haven::read_sav("../data/Lin2023/Study 2.sav") |>
  mutate(Gender = as.character(ifelse(Sex == "男", "Male", ifelse(Sex == "女", "Female", "Other")))) |>
  rename(
    Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
    Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
    Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
    Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
    Affective_touch = Touch, Itch = ITCH)
  • VonMohr (2023) - Study 3
    • https://osf.io/7p9u5/
Code
df9 <- read.csv("../data/VonMohr2023/DataSet_study3.csv")
df9 <- df9[complete.cases(select(df9, starts_with("IAS_"))), ]
df9 <- df9 |>
  mutate(Gender = as.character(ifelse(GENDER == "Man", "Male", ifelse(GENDER == "Woman", "Female", "Other")))) |>
  rename(
    Age=AGE, Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )
  • Makowski
    • Note: First sample has some missing items (No Taste, Cough, Blood_Sugar)
    • Note: Analog scales
    • note df10c = primals data
Code
df10a <- read.csv("https://raw.githubusercontent.com/RealityBending/IllusionGameReliability/main/data/preprocessed_questionnaires.csv") |>
  rename(
    Gender = Sex, Heart = Item_IAS_Interoception_1, Hungry = Item_IAS_Interoception_2,
    Breathing = Item_IAS_Interoception_3, Thirsty = Item_IAS_Interoception_4,
    Temp = Item_IAS_Interoception_5, Sex_arousal = Item_IAS_Interoception_6,
    Urinate = Item_IAS_Elimination_1, Defecate = Item_IAS_Elimination_2,
    Vomit = Item_IAS_Elimination_3, Wind = Item_IAS_Expulsion_1,
    Burp = Item_IAS_Expulsion_2, Sneeze = Item_IAS_Expulsion_3,
    Muscles = Item_IAS_Nociception_1, Bruise = Item_IAS_Nociception_2,
    Pain = Item_IAS_Nociception_3, Affective_touch = Item_IAS_Skin_1,
    Tickle = Item_IAS_Skin_2, Itch = Item_IAS_Skin_3
  ) |>
  filter(!is.na(Urinate))

df10b <- read.csv("https://raw.githubusercontent.com/DominiqueMakowski/PHQ4R/main/study2/data/data.csv") |>
  rename(
    Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21,
    typeSample = Sample
  )

df10c <- read.csv("https://raw.githubusercontent.com/RealityBending/InteroceptionPrimals/refs/heads/main/data/data_participants.csv") |>
  select("Participant" = "participant_id", "Gender", "Age", "Ethnicity", starts_with(c("IAS","MAIA","PI"))) |>
   rename(
    Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  ) |>
  filter(!is.na(Heart)) # participant is outlier and had total IAS scores below 0.4 

df10d <- read.csv("https://raw.githubusercontent.com/RealityBending/InteroceptionScale/refs/heads/main/study2/data/rawdata_participants.csv") |>
  select("Participant", "Gender", "Age", "Ethnicity", starts_with(c("IAS","MAIA","PI", "PHQ4"))) |>
   rename(
    Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  ) |> 
  filter(Participant %in% read.csv("https://raw.githubusercontent.com/RealityBending/InteroceptionScale/refs/heads/main/study2/data/data_participants.csv")$Participant)
Code
# pi_vars <- c("PI_Enticing", "PI_Alive", "PI_Safe", "PI_Good", "PI_Changing", 
#              "PI_Hierarchical", "PI_Understandable", "PI_AboutMe", "PI_Abundant", 
#              "PI_Acceptable", "PI_Beautiful", "PI_Cooperative", "PI_Funny", 
#              "PI_Harmless", "PI_Improvable", "PI_Intentional", "PI_Interconnected", 
#              "PI_Interesting", "PI_Just", "PI_Meaningful", "PI_NeedsMe", 
#              "PI_Pleasurable", "PI_Progressing", "PI_Regenerative", "PI_Stable", 
#              "PI_WorthExploring")
# 
# 
# correlation::correlation(
#   df10c[, "IAS_Total", drop = FALSE],  
#   data2 = select(df10c, all_of(pi_vars)), p_adjust = "none"
# ) |> 
#   # Formatting correlation results
#   mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE))) |>
#   
#   # Releveling factors for visualization
#   mutate(Parameter2 = fct_relevel(Parameter2, rev(pi_vars)),
#          Parameter1 = fct_relevel(Parameter1, "IAS_Total")) |>
#   
#   # Plotting the correlation matrix
#   ggplot(aes(x = Parameter1, y = Parameter2)) +
#   geom_tile(aes(fill = r), color = "white") +  
#   geom_text(aes(label = val), size = 3) +  
#   labs(title = "Correlation Matrix") +
#   scale_fill_gradient2(low = "blue", high = "red", mid = "white", limit = c(-1, 1), midpoint = 0, guide = guide_colourbar(ticks = FALSE)) +  
#   theme_minimal() +
#   theme(
#     legend.title = element_blank(),
#     axis.title.x = element_blank(),
#     axis.title.y = element_blank(),
#     axis.text.x = element_text(hjust = 1)
#   )
  • Poreiro et al.,
    • Study 1: https://osf.io/49wbv
    • Study 2: unpublished
Code
df11a <- haven::read_sav("../data/Giulia/Interoceptive Attention ESM.sav") |> 
  select(AGE, GENDER, contains("IAS_ACC"), -IAS_ACC_20, starts_with(c("MAIA", "TAS20", "DEP", "ANX")),  -ANX_21, -TAS20_21) |> #anx_21 , tas20_21 are attention checks
  mutate_all(as.numeric) |>
  rename(
    Age = AGE, Gender = GENDER,
    Heart = IAS_ACC_1, Hungry = IAS_ACC_2, Breathing = IAS_ACC_3, Thirsty = IAS_ACC_4,
    Urinate = IAS_ACC_5, Defecate = IAS_ACC_6, Taste = IAS_ACC_7, Vomit = IAS_ACC_8, Sneeze = IAS_ACC_9,
    Cough = IAS_ACC_10, Temp = IAS_ACC_11, Sex_arousal = IAS_ACC_12, Wind = IAS_ACC_13, Burp = IAS_ACC_14,
    Muscles = IAS_ACC_15, Bruise = IAS_ACC_16, Pain = IAS_ACC_17, Blood_Sugar = IAS_ACC_18,
    Affective_touch = IAS_ACC_19, Tickle =IAS_ACC_21, Itch = IAS_ACC_22
  ) |>
   mutate(Gender = case_when(
    Gender == 1 ~ "Male",
    Gender == 2 ~ "Female",
    Gender %in% c(3, 5) ~ "Other",
    Gender == 4 ~ "NA"
  )) |> ## compute dimensions for the TAS
  mutate(TAS_DIF = TAS20_1 + TAS20_3 + TAS20_6 + TAS20_7 + TAS20_9 + TAS20_13 + TAS20_14,
         TAS_DDF = TAS20_2 + (1-TAS20_4) + TAS20_11 + TAS20_12 + TAS20_17,
         TAS_EOT = (1-TAS20_5) + TAS20_8 + (1-TAS20_10) + TAS20_15 + TAS20_16 + (1-TAS20_18) + (1-TAS20_19) + TAS20_20)


df11b <- haven::read_sav("../data/Giulia/Sleep and Exteroceptive Interoceptive Sensitivity.sav") |> # 165 participants
  filter(Total_Attention_Fail == 0) |> # 32 participants removed based on failing one or more checks
  select(AGE, GENDER, contains("IAS_ACC"), ) |> 
  mutate_all(as.numeric) |> 
  rename( Age = AGE, Gender = GENDER, 
          Heart = IAS_ACC_1, Hungry = IAS_ACC_2, Breathing = IAS_ACC_3, Thirsty = IAS_ACC_4, Urinate = IAS_ACC_5, Defecate = IAS_ACC_6, Taste = IAS_ACC_7,
          Vomit = IAS_ACC_8, Sneeze = IAS_ACC_9, Cough = IAS_ACC_10, Temp = IAS_ACC_11, Sex_arousal = IAS_ACC_12, Wind = IAS_ACC_13, Burp = IAS_ACC_14,
          Muscles = IAS_ACC_15, Bruise = IAS_ACC_16, Pain = IAS_ACC_17, Blood_Sugar = IAS_ACC_18, Affective_touch = IAS_ACC_19, Tickle =IAS_ACC_20, 
          Itch = IAS_ACC_21 )  |> 
  mutate(Gender = case_when(
    Gender == 1 ~ "Male", 
    Gender == 2 ~"Female", 
    Gender %in% c(3, 5) ~ "Other", 
    Gender == 4 ~ "NA", )) |>
  na.omit()


# Save all
data <- list(df1a=df1a, df1b=df1b, df2=df2, df3=df3, df4=df4, df5=df5, df6=df6, df7a=df7a, df7b=df7b, df7c=df7c, df8a=df8a, df8b=df8b, df9=df9, df10a=df10a, df10b=df10b, df10c=df10c, df11a=df11a, df11b=df11b)
save(data, file = "../data/data.RData")

Participants

  • Sample 1a: Data from Murphy’s (2020) study 1, downloaded from OSF, included 451 participants (Mean age = 25.8, SD = 8.4, range: [18, 69]; Gender: 69.4% women, 29.5% men, 1.11% non-binary).
  • Sample 1b: Data from Murphy’s (2020) study 6, downloaded from OSF, included 375 participants (Mean age = 35.3, SD = 16.9, range: [18, 91]; Gender: 70.1% women, 28.5% men, 1.33% non-binary).
  • Sample 2: Data from Gaggero’s(2020) study, downloaded from OSF, included 814 participants (Mean age = 24.9, SD = 5.3, range: [18, 58], 0.2% missing; Gender: 60.3% women, 39.4% men, 0.25% non-binary).
  • Sample 3: Data from Campos’s(2022) study, downloaded from OSF, included 515 participants (Mean age = 30.7, SD = 10.5, range: [18, 72]; Sex: 0.0% females, 0.0% males, 100.0% other; Gender: 59.6% women, 40.4% men, 0.00% non-binary).
  • Sample 4: Data from Todd’s(2022) study, downloaded from OSF, included 802 participants ().
  • Sample 5: Data from Arslanova (2022) study, downloaded from OSF, included 143 participants (Mean age = 28.5, SD = 7.6, range: [18, 73]; Gender: 45.5% women, 54.5% men, 0.00% non-binary).
  • Sample 6: Data from Brand’s(2022) study, downloaded from OSF, included 619 participants (Mean age = 43.9, SD = 14.5, range: [18, 78]; Gender: 78.7% women, 20.2% men, 1.13% non-binary).
  • Sample 7a: Data from Brand’s(2023) study, downloaded from OSF, included 522 participants (Mean age = 23.4, SD = 6.7, range: [18, 79]; Gender: 79.5% women, 19.7% men, 0.77% non-binary).
  • Sample 7b: Data from Brand’s(2023) study, downloaded from OSF, included 1993 participants (Mean age = 32.0, SD = 12.6, range: [16, 81]; Gender: 77.7% women, 21.7% men, 0.60% non-binary).
  • Sample 7c: Data from Brand’s(2023) study, downloaded from OSF, included 808 participants (Mean age = 27.3, SD = 9.3, range: [18, 72], 0.5% missing; Gender: 68.9% women, 30.2% men, 0.87% non-binary).
  • Sample 8a: Data from Lin’s(2023) study, downloaded from OSF, included 1166 participants (Mean age = 32.5, SD = 8.4, range: [16, 60]; Gender: 57.0% women, 43.0% men, 0.00% non-binary).
  • Sample 8b: Data from Lin’s(2023) study, downloaded from OSF, included 500 participants (Mean age = 37.4, SD = 7.4, range: [20, 60]; Sex: 0.0% females, 0.0% males, 100.0% other; Gender: 56.2% women, 43.8% men, 0.00% non-binary).
  • Sample 9: Data from VonMohr’s(2023) study 3, downloaded from OSF, included 21843 participants (Mean age = 56.5, SD = 14.4, range: [18, 93], 0.2% missing; Gender: 73.2% women, 25.1% men, 1.55% non-binary, 0.15% missing).
  • Sample 10a: Data from [Makowski’s(2023)] study , downloaded from OSF, included 485 participants (Mean age = 30.1, SD = 10.1, range: [18, 73]; Gender: 50.3% women, 49.7% men, 0.00% non-binary; Education: Bachelor, 45.15%; Doctorate, 1.86%; High school, 34.43%; Master, 15.88%; Other, 2.47%; Prefer not to say, 0.21%).
  • Sample 10b: Data from [Makowski’s(2023)] study , downloaded from OSF, included 836 participants (Mean age = 25.1, SD = 11.3, range: [18, 76], 0.1% missing; Gender: 73.8% women, 22.6% men, 2.87% non-binary, 0.72% missing; Education: Bachelor, 22.85%; Doctorate, 2.15%; High School, 63.52%; Master, 6.22%; missing, 0.24%; Other, 5.02%).
  • Sample 10c: Data from [Makowski’s(2023)] study , downloaded from OSF, included 146 participants (Mean age = 21.1, SD = 4.3, range: [18, 50], 0.7% missing; Gender: 80.8% women, 15.1% men, 2.74% non-binary, 1.37% missing).
  • Sample 11a: Data from Poreiro’s(2024) study , included 107 participants (Mean age = 26.8, SD = 9.2, range: [18, 57]; Gender: 74.8% women, 23.4% men, 1.87% non-binary)
  • Sample 11b: Data from [Poreiro’s] study , included 131 participants (Mean age = 30.9, SD = 12.0, range: [18, 60]; Gender: 76.3% women, 22.9% men, 0.76% non-binary)

Total N = 32256.

Code
desc <- function(df){
  summary <- df |> summarise(
      n = n(),
      mean = mean(Age, na.rm = TRUE),
      sd = sd(Age, na.rm = TRUE),
      min = min(Age, na.rm = TRUE),
      max = max(Age, na.rm = TRUE),
      female = sum(Gender == "Female", na.rm = TRUE)
    )
    return(summary)
}

desc_df1a <- desc(df1a)
desc_df1b <- desc(df1b)
desc_df2 <- desc(df2)
desc_df3 <- desc(df3)
#desc_df4 <- desc(df4) # no demographic data available
desc_df5 <- desc(df5_attention)
desc_df6 <- desc(df6)
desc_df7a <- desc(df7a)
desc_df7b <- desc(df7b)
desc_df7c <- desc(df7c)
desc_df8a <- desc(df8a)
desc_df8b <- desc(df8b)
desc_df9 <- desc(df9)
desc_df10a <- desc(df10a)
desc_df10b <- desc(df10b)
desc_df10c <- desc(df10c)
desc_df11a <- desc(df11a)
desc_df11b <- desc(df11b)

desc_total <- rbind(desc_df1a, desc_df1b, desc_df2, desc_df3, desc_df5, desc_df6, desc_df7a, desc_df7b, desc_df7c, desc_df8a, desc_df8b, desc_df9, desc_df10a, desc_df10b, desc_df10c, desc_df11a, desc_df11b)

## calculate weighted mean

total_mean = sum(desc_total$mean * desc_total$n) / sum(desc_total$n)
total_sd = sum(desc_total$sd * desc_total$n) / sum(desc_total$n)
perce_female = sum(desc_total$female)/sum(desc_total$n) *100
Code
library(gt)

# APA style ####
gt_apastyle <- function(gt_table, font.size=12) {
  gt_table  |> 
    gt::opt_table_lines(extent = "none") %>%
    gt::tab_options(
      heading.border.bottom.width = 2,
      heading.border.bottom.color = "black",
      heading.border.bottom.style = "solid",
      table.border.top.color = "black",
      table.border.top.style = "solid",
      table.border.top.width = 2,  
      table_body.hlines.color = "white",
      table_body.border.top.color = "black",
      table_body.border.top.style = "solid",
      table_body.border.top.width = 2,
      heading.title.font.size = font.size,
      table.font.size = font.size,
      heading.subtitle.font.size = font.size,
      table_body.border.bottom.color = "black",
      table_body.border.bottom.width = 2,
      table_body.border.bottom.style = "solid",
      column_labels.border.bottom.color = "black",
      column_labels.border.bottom.style = "solid",
      column_labels.border.bottom.width = 1,
      latex.use_longtable = FALSE
    ) |> 
      gt::opt_table_font(font = "times")
}

make_age <- function(age) {
  age <- as.numeric(age) 
  
  mean(age, na.rm = TRUE) |> 
    insight::format_value(digits=1) |> 
    paste0(" ± ", 
           insight::format_value(sd(age, na.rm = TRUE), digits=1)) 
} 

# Create the dataframe
table <- data.frame(
  Reference = c('Murphy et al., (2020)', '', '', 'Gaggero et al., (2021)', 'Campos et al., (2022)', 'Todd et al., (2022)', 'Arslanova et al., (2022)', 'Brand et al., (2022)', 'Brand et al., (2023)', '', '', '', 'Lin et al., (2023)', '', '', 'VonMohr et al., (2023)', 'Makowski et al., (2023a)', 'Makowski et al., (2023b)', 'Makowski et al., (2023c)', 'Poerio et al., (2024)', 'Poerio et al., unpublished', 'Total'),
  Sample = c('', 'Sample 1a', 'Sample 1b', 'Sample 2', 'Sample 3', 'Sample 4', 'Sample 5', 'Sample 6', '', 'Sample 7a', 'Sample 7b', 'Sample 7c', '', 'Sample 8a', 'Sample 8b', 'Sample 9', 'Sample 10', 'Sample 11', 'Sample 12', 'Sample 13', 'Sample 14', 'Sample 15'),   
  Language = c('', 'English', 'English' , 'English and Italian', 'Portuguese', 'English', 'English', 'German', '', 'German', 'German', 'German', '', 'Chinese', 'Chinese', 'English', 'English', 'English', 'English', 'English', 'English',''),
  N = c('', nrow(df1a), nrow(df1b), nrow(df2), nrow(df3), nrow(df4), nrow(df5_attention), nrow(df6),'', nrow(df7a), nrow(df7b), 802,'', nrow(df8a), nrow(df8b), nrow(df9), nrow(df10a), nrow(df10b), nrow(df10c), nrow(df11a), nrow(df11b), 32214),
  'Difference' = c('','', '', '', '', '', '', '', '', '', '', '', '', 'Collapsed "Itch" and "Tickling"', 'Collapsed "Itch" and "Tickling"', '', 'Analog scales', 'Analog scales', 'Analog scales', '', '', ''),
  Age = c('', make_age(df1a$Age), make_age(df1b$Age), make_age(df2$Age), make_age(df3$Age), "48.6.6 ± 14.1*", make_age(df5_attention$Age), make_age(df6$Age),'', make_age(df7a$Age), make_age(df7b$Age), make_age(df7c$Age),'', make_age(df8a$Age), make_age(df8b$Age), make_age(df9$Age), make_age(df10a$Age), make_age(df10b$Age), make_age(df10c$Age), make_age(df11a$Age), make_age(df11b$Age), '48.6 ± 13.1'),
  Range = c('', '18-69', '18-91', '18-58', '18-72', '18-92*', '18-73', '18-78', '', '18-79', '16-81', '18-72','', '16-60', '20-60', '18-93', '18-73', '17-76','18-50', '18-57', '18-60', '17-93'),
  Female_Percentage = c('', '69.4%', '70.1%', '60.3%', '59.6%', '50%*', '46.8%', '78.7%', '', '79.5%', '77.7%', '68.9%', '', '57.0%', '56.2%', '73.2%', '50.3%', '53.0%', '76%', '74.8%', '75.9%', '71.6%'),
  Availability= c('osf.io/3m5nh', '', '', 'osf.io/5x9sg', 'osf.io/j6ef3', 'osf.io/ms354', 'osf.io/mp3cy', 'osf.io/xwz6g', 'osf.io/3f2h6', '', '', '','osf.io/3eztd','', '', 'osf.io/7p9u5', 'github.com/RealityBending/IllusionGameReliability', 'github.com/DominiqueMakowski/PHQ4R', 'github.com/RealityBending/InteroceptionPrimals', 'osf.io/49wbv', '','')
) 

table_apa <- table |> 
  gt() |>
  cols_align(align = c("right"), columns = "Age") |> 
  cols_label(Age = "Age (Mean  ± SD)", Female_Percentage = "Female %") |> 
  tab_footnote("* Information taken from the sample description of relevant paper rather than recomputed.")
  
gt_apastyle(table_apa, font.size=12)
gtsave(gt_apastyle(table_apa, font.size=9), "figures/table1.tex")

# gtsave(table_apa, "figures/table1.png")

Item Selection

Distribution

The items with the differing distribution pattern (with a lower mode around 2/5) are “Affective Touch”, “Blood Sugar” and “Bruise”.

Code
vars <- c(
  "Heart", "Hungry", "Breathing", "Thirsty", "Urinate", "Defecate", "Taste", "Vomit", "Sneeze", "Cough", "Temp",
  "Sex_arousal", "Wind", "Burp", "Muscles", "Bruise", "Pain", "Blood_Sugar", "Affective_touch", "Tickle", "Itch"
)

dens1a <- estimate_density(select(df1a, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 1a")
dens1b <- estimate_density(select(df1b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 1b")
dens2 <- estimate_density(select(df2, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 2")
dens3 <- estimate_density(select(df3, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 3")
dens4 <- estimate_density(select(df4, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 4")
dens5 <- estimate_density(select(df5, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 5")
dens6 <- estimate_density(select(df6, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 6")
dens7a <- estimate_density(select(df7a, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7a")
dens7b <- estimate_density(select(df7b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7b")
dens7c <- estimate_density(select(df7c, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7c")
dens8a <- estimate_density(select(df8a, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
  mutate(Sample = "Sample 8a")
dens8b <- estimate_density(select(df8b, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
  mutate(Sample = "Sample 8b")
dens9 <- estimate_density(select(df9, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 9")
dens10a <- estimate_density(select(df10a, all_of(setdiff(vars, c("Taste", "Cough", "Blood_Sugar")))), method = "kernSmooth") |>
  mutate(Sample = "Sample 10a",
         x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens10b <- estimate_density(select(df10b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 10b",
         x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens10c <- estimate_density(select(df10c, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 10c",
         x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens11a <- estimate_density(select(df11a, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 11a")
dens11b <- estimate_density(select(df11b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 11b")


p1_data <- rbind(dens1a, dens1b, dens2, dens3, dens4, dens5, dens6, dens7a, dens7b, dens7c, dens8a, dens8b, dens9, dens10a, dens10b, dens10c, dens11a, dens11b) 

linetype <- setNames(rep("solid", length(vars)), vars)
linetype["Affective touch"] <- "dotted"
linetype["Blood sugar"] <- "dashed"
linetype["Sex arousal"] <- "solid"
linetype["Bruise"] <- "dashed"

# proportion of non-1 and non-5 values 
p1_data |>
  filter(x > 1 & x < 5) |>  # Keep only x values strictly between 1 and 5
  summarise(proportion = sum(y) / sum(p1_data$y))  # Normalize by total density
  proportion
1  0.9994812
Code
p1 <- p1_data |>
  mutate(Sample = fct_relevel(Sample, "Sample 1a", "Sample 1b", "Sample 2", "Sample 3", "Sample 4", "Sample 5", "Sample 6", "Sample 7a", "Sample 7b", "Sample 7c", "Sample 8a", "Sample 8b", "Sample 9", "Sample 10a", "Sample 10b", "Sample 10c", "Sample 11a", "Sample 11b"),
         Parameter = ifelse(Parameter == "Affective_touch", "Affective touch", Parameter),
         Parameter = ifelse(Parameter == "Blood_Sugar", "Blood sugar", Parameter),
         Parameter = ifelse(Parameter == "Sex_arousal", "Sex arousal", Parameter)) |>
  # mutate(Dashed = ifelse(Parameter %in% c("Affective_touch", "Blood_Sugar", "Bruise"), TRUE, FALSE)) |> 
  ggplot(aes(x = x, y = y, color = Parameter)) +
  geom_line(aes(linetype = Parameter), linewidth = 0.7) +
  scale_color_material_d() +
  scale_linetype_manual(values = linetype)  +
  labs(x = "Response", y = "Distribution", color = "Item", linetype = "Item", title = "Item Distribution") +
  guides(color = guide_legend(ncol = 1)) +
  facet_wrap(~Sample, scales = "free_y", nrow = 3) +
  theme_minimal() +
  theme(axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        plot.title = element_text(face="bold"))

p1

Code
data1a <- normalize(select(df1a, all_of(dens1a$Parameter)), range = c(1, 5))
data1b <- normalize(select(df1b, all_of(dens1b$Parameter)), range = c(1, 5))
data2 <- normalize(select(df2, all_of(dens2$Parameter)), range = c(1, 5))
data3 <- normalize(select(df3, all_of(dens3$Parameter)), range = c(1, 5))
data4 <- normalize(select(df4, all_of(dens4$Parameter)), range = c(1, 5))
data5 <- normalize(select(df5, all_of(dens5$Parameter)), range = c(1, 5))
data6 <- normalize(select(df6, all_of(dens6$Parameter)), range = c(1, 5))
data7a <- normalize(select(df7a, all_of(dens7a$Parameter)), range = c(1, 5))
data7b <- normalize(select(df7b, all_of(dens7b$Parameter)), range = c(1, 5))
data7c <- normalize(select(df7c, all_of(dens7c$Parameter)), range = c(1, 5))
data8a <- normalize(select(df8a, all_of(dens8a$Parameter)), range = c(1, 5))
data8b <- normalize(select(df8b, all_of(dens8b$Parameter)), range = c(1, 5))
data9 <- normalize(select(df9, all_of(dens9$Parameter)), range = c(1, 5))
data10a <- select(df10a, all_of(dens10a$Parameter))
data10b <- select(df10b, all_of(dens10b$Parameter))
data10c <- select(df10c, all_of(dens10c$Parameter))
data11a <- normalize(select(df11a, all_of(dens11a$Parameter)), range = c(1, 5))
data11b <- na.omit(normalize(select(df11b, all_of(dens11b$Parameter)), range = c(1, 5))) # remove NAS
row.names(data11b) <- NULL



data_all <- rbind(
  data1a, data1b, data2, data3, data4, data5, data6, data7a, data7b, data7c,
  mutate(data8a, Tickle = NA), 
  mutate(data8b, Tickle = NA), data9,
  mutate(data10a, Taste = NA, Cough = NA, Blood_Sugar = NA),  
  data10b,  data10c, data11a, data11b
) 

Correlations

An overall positive intercorrelation pattern, with no clear structure emerging.

Code
make_correlation <- function(df) {
  correlation::correlation(df, redundant = TRUE) |>
    correlation::cor_sort() |>
    # correlation::cor_lower() |>
    mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE)),
           r = ifelse(Parameter1 == Parameter2, NA, r),
           lbl = ifelse(Parameter1 == Parameter2, "", format_value(r, zap_small = TRUE, lead_zero = FALSE)),
           lbl = ifelse(p > .05, "", lbl),
           Param1 = fct_relevel(str_replace(Parameter1, "_", " "), str_replace(levels(Parameter1), "_", " ")),
           Param2 = fct_relevel(str_replace(Parameter2, "_", " "), str_replace(levels(Parameter2), "_", " "))) |>
    # mutate(Parameter2 = fct_rev(Parameter2)) |>
    ggplot(aes(x = Param1, y = Param2)) +
    geom_tile(aes(fill = r), color = "white") +
    geom_text(aes(label = lbl), size = 3) +
    labs(title = "Correlation Matrix", subtitle = paste0("N = ", nrow(df))) +
    scale_fill_gradientn(
      colors = c("white", "#FFF9C4", "#FFF59D", "#FFEB3B", "#FFCA28", "#FF9800", "#FF5722", "#F44336", "#E91E63", "#9C27B0", "#673AB7", "#512DA8", "#311B92"),
      na.value = "#311B92",
      limits = c(0, 0.8),
      guide = guide_colorbar(ticks.colour = NA)
    ) +
    # scale_fill_metro_c(limit = c(0, 0.75), guide = guide_colourbar(ticks = FALSE)) +
    theme_minimal() +
    theme(
      legend.title = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.x = element_text(angle = 45, hjust = 1),
      plot.title = element_text(face = "bold")
    )
}

p2 <- make_correlation(data_all)
p2

Code
fig1 <- p1 / p2
ggsave("figures/Figure1.png", fig1, width=10, height=12, dpi=300, bg="white")
Code
make_correlation(data1a)

Code
make_correlation(data1b)

Code
make_correlation(data2)

Code
make_correlation(data3)

Code
make_correlation(data4)

Code
make_correlation(data5)

Code
make_correlation(data6)

Code
make_correlation(data7a)

Code
make_correlation(data7b)

Code
make_correlation(data7c)

Code
make_correlation(data8a)

Code
make_correlation(data8b)

Code
make_correlation(data9)

Code
make_correlation(data10a)

Code
make_correlation(data10b)

Code
make_correlation(data10c)

Code
make_correlation(data11a)

Code
make_correlation(data11b)

Unique Variable Analysis (UVA)

Unique Variable Analysis (Christensen, Garrido, & Golino, 2023) uses the weighted topological overlap measure (Nowick et al., 2009) on an estimated network. Values greater than 0.25 are determined to have considerable local dependence (i.e., redundancy) that should be handled (variables with the highest maximum weighted topological overlap to all other variables (other than the one it is redundant with) should be removed).

Code
uva0 <- EGAnet::UVA(data = data_all, cut.off = 0.3)
uva0
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.363

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.290
 Urinate Defecate 0.265

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
 Sneeze     Cough 0.230
  Heart Breathing 0.225
 Hungry   Thirsty 0.217
Code
uva0$keep_remove
$keep
[1] "Itch"

$remove
[1] "Tickle"
Code
uva1a <- EGAnet::UVA(data = data1a, cut.off = 0.3)
uva1a
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.286
   Wind   Burp 0.275

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Sneeze    Cough 0.244
 Urinate Defecate 0.241
Code
uva1a$keep_remove
NULL
Code
uva1b <- EGAnet::UVA(data = data1b, cut.off = 0.3)
uva1b
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.350
 Sneeze  Cough 0.317
   Wind   Burp 0.309

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva1b$keep_remove
$keep
[1] "Cough" "Burp"  "Itch" 

$remove
[1] "Sneeze" "Wind"   "Tickle"
Code
uva2 <- EGAnet::UVA(data = data2, cut.off = 0.3)
uva2
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
 Tickle      Itch 0.283
  Heart Breathing 0.252

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Sneeze    Cough 0.230
 Urinate Defecate 0.223
  Hungry  Thirsty 0.200
Code
uva2$keep_remove
NULL
Code
uva3 <- EGAnet::UVA(data = data3, cut.off = 0.3)
uva3
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.445
 Sneeze  Cough 0.318

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
   Wind   Burp 0.256

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i      node_j   wto
  Bruise Blood_Sugar 0.219
 Urinate    Defecate 0.217
   Heart   Breathing 0.208
Code
uva3$keep_remove
$keep
[1] "Sneeze" "Itch"  

$remove
[1] "Cough"  "Tickle"
Code
uva4 <- EGAnet::UVA(data = data4, cut.off = 0.3)
uva4
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.415
 Tickle      Itch 0.351

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
  Hungry  Thirsty 0.293
 Urinate Defecate 0.282
    Wind     Burp 0.275
  Sneeze    Cough 0.251

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i node_j   wto
 Muscles   Pain 0.205
Code
uva4$keep_remove
$keep
[1] "Heart" "Itch" 

$remove
[1] "Breathing" "Tickle"   
Code
uva5 <- EGAnet::UVA(data = data5, cut.off = 0.3)
Warning: Some variables did not belong to a dimension: Blood_Sugar, Affective_touch 

Use caution: These variables have been removed from the TEFI calculation
Code
uva5
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

    node_i node_j   wto
    Tickle   Itch 0.247
    Bruise   Itch 0.246
 Breathing  Vomit 0.207
Code
uva5$keep_remove
NULL
Code
uva6 <- EGAnet::UVA(data = data6, cut.off = 0.3)
uva6
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.416
  Sneeze    Cough 0.332
    Wind     Burp 0.314
  Tickle     Itch 0.303

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i  node_j   wto
 Hungry Thirsty 0.274

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.237
Code
uva6$keep_remove
$keep
[1] "Urinate" "Sneeze"  "Wind"    "Itch"   

$remove
[1] "Defecate" "Cough"    "Burp"     "Tickle"  
Code
uva7a <- EGAnet::UVA(data = data7a, cut.off = 0.3)
uva7a
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Sneeze  Cough 0.434
 Tickle   Itch 0.376

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.295
 Urinate Defecate 0.294

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva7a$keep_remove
$keep
[1] "Sneeze" "Tickle"

$remove
[1] "Cough" "Itch" 
Code
uva7b <- EGAnet::UVA(data = data7b, cut.off = 0.3)
uva7b
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.317
  Tickle     Itch 0.311

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.281
 Sneeze     Cough 0.261

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.237
Code
uva7b$keep_remove
$keep
[1] "Defecate" "Itch"    

$remove
[1] "Urinate" "Tickle" 
Code
uva7c <- EGAnet::UVA(data = data7c, cut.off = 0.3)
uva7c
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.457
 Sneeze  Cough 0.302

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.280
 Urinate Defecate 0.254

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
 Hungry   Thirsty 0.242
  Vomit    Sneeze 0.224
  Heart Breathing 0.218
Code
uva7c$keep_remove
$keep
[1] "Cough" "Itch" 

$remove
[1] "Sneeze" "Tickle"
Code
uva8a <- EGAnet::UVA(data = data8a, cut.off = 0.3)
uva8a
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.237
   Heart Breathing 0.218
  Hungry   Thirsty 0.213
Code
uva8a$keep_remove
NULL
Code
uva8b <- EGAnet::UVA(data = data8b, cut.off = 0.3)
uva8b
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.277
   Heart Breathing 0.267

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.206
Code
uva8b$keep_remove
NULL
Code
uva9 <- EGAnet::UVA(data = data9, cut.off = 0.3)
uva9
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
  Tickle     Itch 0.379
    Wind     Burp 0.373
 Urinate Defecate 0.341

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.285
 Sneeze     Cough 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i  node_j  wto
 Hungry Thirsty 0.24
Code
uva9$keep_remove
$keep
[1] "Defecate" "Wind"     "Itch"    

$remove
[1] "Urinate" "Burp"    "Tickle" 
Code
uva10a <- EGAnet::UVA(data = data10a, cut.off = 0.3)
uva10a
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.297

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.209
Code
uva10a$keep_remove
NULL
Code
uva10b <- EGAnet::UVA(data = data10b, cut.off = 0.3)
uva10b
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva10b$keep_remove
NULL
Code
uva10c <- EGAnet::UVA(data = data10c, cut.off = 0.3)
Warning: Some variables did not belong to a dimension: Cough 

Use caution: These variables have been removed from the TEFI calculation
Code
uva10c
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva10c$keep_remove
NULL
Code
uva11a <- EGAnet::UVA(data = data11a, cut.off = 0.3)
uva11a
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.330
   Heart Breathing 0.314
    Wind      Burp 0.302

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i  node_j   wto
 Muscles    Pain 0.279
  Tickle    Itch 0.277
  Hungry Thirsty 0.255

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
  Vomit Sneeze 0.201
Code
uva11a$keep_remove
$keep
[1] "Heart"    "Defecate" "Burp"    

$remove
[1] "Breathing" "Urinate"   "Wind"     
Code
uva11b <- EGAnet::UVA(data = data11b, cut.off = 0.3)
uva11b
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
 Hungry Breathing 0.252

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i  node_j   wto
   Wind    Burp 0.214
   Temp Muscles 0.202
Code
uva11b$keep_remove
NULL
Code
rez_uva <- c(uva1a$keep_remove$remove,
  uva1b$keep_remove$remove,
  uva2$keep_remove$remove,
  uva3$keep_remove$remove,
  uva4$keep_remove$remove,
  uva5$keep_remove$remove,
  uva6$keep_remove$remove,
  uva7a$keep_remove$remove,
  uva7b$keep_remove$remove,
  uva7c$keep_remove$remove,
  uva8a$keep_remove$remove,
  uva8b$keep_remove$remove,
  uva9$keep_remove$remove,
  uva10a$keep_remove$remove,
  uva10b$keep_remove$remove,
  uva10c$keep_remove$remove,
  uva11a$keep_remove$remove,
  uva11b$keep_remove$remove)

sort(table(rez_uva), decreasing = TRUE) / length(rez_uva) 
rez_uva
    Tickle      Cough    Urinate  Breathing       Burp     Sneeze       Wind 
0.30434783 0.13043478 0.13043478 0.08695652 0.08695652 0.08695652 0.08695652 
  Defecate       Itch 
0.04347826 0.04347826 

Structure Refinement

Discarded the following items: - Tickle: not present in the same dataset and consistently flagged as redundant in UVA analysis.

Code
data_all <- select(data_all, -Tickle)
data1a <- select(data1a, -Tickle)
data1b <- select(data1b, -Tickle)
data2 <- select(data2, -Tickle)
data3 <- select(data3, -Tickle)
data4 <- select(data4, -Tickle)
data5 <- select(data5, -Tickle)
data6 <- select(data6, -Tickle)
data7a <- select(data7a, -Tickle)
data7b <- select(data7b, -Tickle)
data7c <- select(data7c, -Tickle)
# data8a <- select(data8a, -Tickle)
# data8b <- select(data8b, -Tickle)
data9 <- select(data9, -Tickle)
data10a <- select(data10a, -Tickle)
data10b <- select(data10b, -Tickle)
data10c <- select(data10c, -Tickle)
data11a <- select(data11a, -Tickle)
data11b <- select(data11b, -Tickle)

Initial

Code
colors <- c("Heart" = "#F44336", "Breathing" = "#FF5722",
            "Hungry" = "#FF9800", "Thirsty" = "#FFC107",
            "Burp" = "#4CAF50", "Wind" = "#009688",
            "Urinate" = "#63824C", "Defecate" = "#795548",
            "Cough" = "#8BC34A", "Sneeze" = "#CDDC39", 
            "Bruise" = "#673AB7", "Blood Sugar" = "#3F51B5",
            "Muscles" = "#2196F3", "Pain" = "#00BCD4",
            "Sex arousal" = "#FF4081", "Temp" = "#9C27B0",
            "Vomit" = "#76FF03", "Taste" = "#00E676",
            "Affective touch" = "#FF1744", "Itch" = "#D500F9")
Code
make_hclust <- function(data) {
  rez_pvclust <- pvclust::pvclust(data,
                 method.hclust = "complete",
                 method.dist = "correlation",
                 nboot = 1000, quiet=TRUE, parallel=TRUE)
  
  
  dendrogram <- as.dist(1 - cor(data, use = "pairwise.complete.obs")) |> 
  hclust(method = "complete") |> 
  tidygraph::as_tbl_graph() 


  # Process Nodes
  nodes <- as.list(dendrogram)$nodes |> 
    mutate(
      Size = ifelse(label != "", 10, NA),
      Item = str_replace(label, "_", " "),
      idx = 1:nrow(as.list(dendrogram)$nodes))
  
  # Central node 
  nodes[nodes$height == max(nodes$height), c("Item", "Size")] <- data.frame(Item="Central", Size=15)
  
  # Process Edges
  edges <- as.list(dendrogram)$edges
  edges$linewidth = datawizard::rescale(nodes[edges$from, ]$height, to = c(0.1, 1))
  
  p <- tidygraph::tbl_graph(nodes = nodes, edges = edges) |> 
    ggraph(layout = "dendrogram", circular = TRUE) +
    # geom_edge_diagonal(strength = 0.7, linewidth = 1) +
    geom_edge_elbow2(aes(edge_width=linewidth), color="#212121") +
    geom_node_point(aes(filter=Item %in% c("Central"), size = Size), color="#212121") +
    # geom_node_point(aes(filter=Group %in% c("Visceroception", "Awareness", "Deficit"), color=Group, size = Size)) +
    geom_node_text(aes(
      label = ifelse(Item != "", Item, NA),
      x = x * 1.10,
      y = y * 1.10,
      filter = label != "",
      angle = ifelse(
        x >= 0,
        asin(y) * 360 / 2 / pi,
        360 - asin(y) * 360 / 2 / pi
      ),
      hjust = ifelse(
        x >= 0, 0, 1
      ))) +
    geom_node_point(aes(filter = label != "", color=Item, size=Size), alpha = 1) +
    # geom_node_text(aes(label=idx)) +  # Debug
    scale_edge_width_continuous(range=c(1, 3), guide = "none") +
    scale_size_continuous(range=c(10, 10), guide = "none") +
    scale_color_manual(values = colors,
                       breaks = names(colors)) +
    ggtitle("Hierarchical Clustering Analysis (HCA)", subtitle = "Method = Correlation") +
    coord_equal(clip = "off", xlim = c(-1.25, 1.25), ylim = c(-1.25, 1.25)) +
    theme_void() + 
    # guides(color = guide_legend(override.aes = list(size = c(7.5, 5, 5, 3.5, 7.5, 5, 5, 3.5, 3.5, 3.5, 7.5, 5, 5, 3.5)))) +
    theme(legend.text = ggtext::element_markdown(),
          legend.title = element_blank(),
          legend.position = "none",
          # plot.title = element_blank(), 
          # plot.subtitle = element_blank()
          plot.title = element_text(face="bold"))

  
  list(pvclust = rez_pvclust, p = p)
}


rez_hclust <- make_hclust(data_all)
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
rez_hclust$p

Code
make_ega <- function(data, cols = c("C1" = "#F44336", "C2" = "#FF9800", "C3" = "#795548", "C4" = "#CDDC39", "C5" = "#2196F3", "C6" = "#009688", "C7" = "#9C27B0")) {
  rez <- EGAnet::bootEGA(
    data = data,
    seed = 123,
    model = "glasso",
    algorithm = "leiden",
    EGA.type = "hierEGA",
    type = "resampling",
    plot.itemStability = FALSE,
    verbose = FALSE, 
    allow.singleton = TRUE
  )
  
  table <- EGAnet::net.loads(rez$EGA$lower_order)$std |> 
    as.data.frame() |> 
    rownames_to_column("Item") |>
    gt::gt() |>
    gt::tab_header(title = "EGA Loadings") |>
    gt::data_color(
      columns = -Item,
      method = "numeric",
      colors = scales::col_numeric(
        palette = c("red", "white", "green"),
        domain = c(-1, 0, 1)
      )) |> 
    gt::fmt_auto()
  
    
  nodes <- rez$EGA$lower_order$dim.variables |> 
    rename(name = items) |> 
    mutate(dimension = paste0("C", dimension),
           size = apply(EGAnet::net.loads(rez$EGA$lower_order)$std, 1, max)[name])
  
  loadings <- rez$EGA$lower_order$network |> 
    as.data.frame() |> 
    rownames_to_column("to") |>
    pivot_longer(-c(to), names_to = "from", values_to = "Loading") |> 
    filter(Loading > quantile(Loading, 1/3)) 
  
  g <- tidygraph::tbl_graph(nodes = nodes, edges = loadings, directed = FALSE) |> 
    mutate(name = str_replace(name, "_", " "),
           name = ifelse(name == "Sex arousal", "Sex", name)) 
  
  set.seed(1)
  layout <- ggraph::create_layout(g, layout = "fr", weights = abs(loadings$Loading))
  
  xrange <- max(layout$x) - min(layout$x)
  yrange <- max(layout$y) - min(layout$y)
  xmin <- min(layout$x) - xrange * 0.05
  xmax <- max(layout$x) + xrange * 0.05
  ymin <- min(layout$y) - yrange * 0.05
  ymax <- max(layout$y) + yrange * 0.05
  
  p <- layout |> 
    ggraph() +
    geom_edge_bend(aes(edge_width = Loading, edge_alpha = Loading), color = "#212121",
                   strength = 0.3) +
    geom_node_point(aes(size = size, color = dimension), alpha = 0.95) +
    geom_node_text(aes(label = name), size = 3, color = "white", fontface = "bold") +
    scale_size_continuous(range = c(23, 28)) +
    scale_edge_width(range = c(0.3, 4)) +
    scale_edge_alpha(range = c(0.1, 0.9)) +
    scale_color_manual(values = cols) +
    labs(title = "Exploratory Graph Analysis (EGA)", subtitle = "Method = Leiden") +
    theme_void() +
    theme(legend.position = "none",
          plot.title = element_text(face="bold")) +
    coord_cartesian(xlim =c(xmin, xmax), ylim = c(ymin, ymax)) 
    
  list(rez = rez, table = table, p = p)
}

rez_ega <- make_ega(data_all, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]]))
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.524 0.038 0.009 0.067 0.038 0.019 0.114
Breathing 0.524 0.193 0.046 0.185 0.055 0.002 0    
Hungry 0.133 0.522 0.07  0.124 0     0.003 0.037
Thirsty 0.095 0.522 0.124 0.105 0.009 0.011 0.01 
Urinate 0.021 0.179 0.563 0.15  0.038 0.015 0    
Defecate 0.046 0.057 0.563 0.156 0.113 0.042 0    
Pain 0.037 0.109 0.055 0.387 0.046 0     0.138
Temp 0.026 0.089 0.086 0.365 0.127 0.028 0.007
Muscles 0.103 0.024 0.039 0.359 0.058 0.089 0.127
Sex_arousal 0.073 0.065 0.071 0.298 0.079 0.078 0    
Affective_touch 0.037 0     0.034 0.276 0.05  0     0.207
Taste 0.07  0.03  0.062 0.225 0.116 0.039 0.055
Sneeze 0.03  0     0.051 0.108 0.569 0.081 0.047
Cough 0.038 0.012 0.012 0.136 0.446 0.162 0.086
Vomit 0.047 0     0.094 0.188 0.293 0.043 0.042
Wind 0.008 0.009 0.055 0.116 0.104 0.572 0.059
Burp 0.019 0.01  0.005 0.099 0.184 0.572 0.111
Bruise 0.03  0     0     0.186 0.048 0.081 0.44 
Blood_Sugar 0.067 0.049 0     0.077 0.02  0.025 0.387
Itch 0.019 0     0     0.138 0.076 0.033 0.298
Code
rez_ega$p

Code
make_efa <- function(data, n = 3, cols = c("F1" = "#009688", "F2" = "#795548", "F3" = "red")) {
  n_fac <- n_factors(data)

  rez_efa <- factor_analysis(data, n = n, sort = TRUE, rotation = "oblimin")
  
  # Exctract loadings
  loadings <- attributes(rez_efa)$loadings_long |> 
    select(from = Component, to=Variable, Loading) |> 
    mutate(Type = "Loading",
           to = str_replace(to, "_", " "),
           from = str_replace(from, "MR", "F")) 
  
  # Reorder factors
  fct_order <- c("Pain", "Muscles", "Blood Sugar", "Bruise",
                 "Affective touch", "Sex arousal", "Temp", "Itch", "Tickle", 
                 "Thirsty", "Hungry", "Breathing", "Heart", "Defecate", "Urinate",
                 "Sneeze", "Cough",  "Wind", "Burp", "Vomit", "Taste")
  loadings <- loadings |> 
    mutate(to = fct_relevel(to, fct_order[fct_order %in% unique(loadings$to)])) |> 
    arrange(to) |> 
    mutate(to = as.character(to))

  # Get correlations between factors
  phi <- attributes(rez_efa)$model$Phi
  phi[upper.tri(phi, diag = TRUE)] <- NA
  phi <- phi |> 
    as.data.frame() |> 
    rownames_to_column("from") %>%
    pivot_longer(-from, names_to = "to", values_to = "Loading") |> 
    filter(!is.na(Loading)) |>
    mutate(Type = "Correlation",
           from = str_replace(from, "MR", "F"),
           to = str_replace(to, "MR", "F")) 
  
  # Reverse so that arcs are on the same side
  phi[paste(phi$from, phi$to) %in% c("F3 F1"), c("from", "to")] <- data.frame(from = "F1", to = "F3")
  phi[paste(phi$from, phi$to) %in% c("F2 F1"), c("from", "to")] <- data.frame(from = "F1", to = "F2")
  phi[paste(phi$from, phi$to) %in% c("F3 F2"), c("from", "to")] <- data.frame(from = "F2", to = "F3")
  phi[paste(phi$from, phi$to) %in% c("F4 F1"), c("from", "to")] <- data.frame(from = "F1", to = "F4")
  phi[paste(phi$from, phi$to) %in% c("F4 F2"), c("from", "to")] <- data.frame(from = "F2", to = "F4")
  phi[paste(phi$from, phi$to) %in% c("F4 F3"), c("from", "to")] <- data.frame(from = "F3", to = "F4")
  
  loadings <- rbind(loadings, phi)
  
  # Make layout
  factor_names <- unique(loadings$from[loadings$Type == "Loading"])
  item_names <- unique(loadings$to[loadings$Type == "Loading"])
  
  nodes <- tibble(name = unique(c(loadings$from, loadings$to))) |>
    mutate(type = case_when(
      name %in% factor_names ~ "Factor",
      name %in% item_names ~ "Item",
      TRUE ~ "Other"
    ))
  
  y_position <- seq(0.8, 0.2, length.out = length(factor_names))
  
  layout <- nodes |>
    mutate(x = ifelse(type == "Item", 0, 1),
           y = case_when(
            name == "F1" ~ y_position[1],
            name == "F2" ~ y_position[2],
            name == "F3" ~ y_position[3],
            name == "F4" ~ y_position[4],
            name == "F5" ~ y_position[5],
            name == "F6" ~ y_position[6],
            .default = NA
          ))
  layout[layout$type == "Item", "y"] <- seq(0, 1, length.out = sum(layout$type == "Item"))
  
  
  g <- tidygraph::as_tbl_graph(loadings, directed = TRUE) |>
    left_join(layout, by = c("name"))
  
  p_efa <- ggraph(g, layout = as.data.frame(layout)[c("x", "y")]) + 
    geom_edge_link(
      aes(
        edge_width = abs(Loading),
        edge_alpha = abs(Loading),
        # edge_color = Loading,
        filter = Type == "Loading"
        # label = insight::format_value(Loading, lead_zero = FALSE, zap_small = TRUE)
      ),
      color = "#212121",
      arrow = arrow(length = unit(2.5, "mm"), type = "closed"),
      end_cap = circle(5, 'mm')
    ) +
      geom_edge_arc(
      aes(
        edge_width = abs(Loading),
        filter = Type == "Correlation",
        label = insight::format_value(Loading, lead_zero = FALSE, zap_small = TRUE)
      ),
      edge_color = "forestgreen",
      angle_calc = "along",
      label_dodge = unit(2.5, "mm"),
      force_flip = TRUE,
      strength = 0.2  # controls curvature
    ) +
    geom_node_label(data=filter(layout, type == "Item"), 
                    aes(fill = name), 
                    label = paste(rep(" ", 20), collapse = ""),
                    size = 3, 
                    hjust = 0.5, 
                    nudge_x = -0.05,
                    label.r = unit(0.2, "lines"),
                    label.padding = unit(0.5, "lines"),
                    label.size = 0) +
    geom_node_label(data=filter(layout, type == "Item"), 
                    aes(label = name), 
                    fill = NA,
                    size = 3, 
                    hjust = 0.5, 
                    vjust = 0.5,
                    color = "white",
                    fontface = "bold",
                    nudge_x = -0.05,
                    label.size = 0) +
    geom_node_point(data=filter(layout, type == "Factor"), 
                    aes(color = name), 
                    size = 12) +
    geom_node_text(data=filter(layout, type == "Factor"), 
                   aes(label = name), 
                   size = 4, 
                   color = "white", 
                   fontface = "bold") +
    scale_edge_width(range = c(0.3, 2)) +
    scale_edge_alpha(range = c(0.1, 1)) +
    scale_fill_manual(values = colors) +
    scale_color_manual(values = cols) +
    labs(title = "Exploratory Factor Analysis (EFA)", subtitle = "Method = Oblimin") +
    theme_void() +
    theme(legend.position = "none",
          plot.title = element_text(face="bold")) +
    coord_cartesian(xlim = c(-0.08, 1.08)) 

  list(n_fac=n_fac, efa=rez_efa, p=p_efa)
}


rez_efa <- make_efa(data_all, n = 3, cols = c("F1" = colors[["Wind"]], 
                                              "F2" = colors[["Defecate"]], 
                                              "F3" = colors[["Heart"]]))
Loading required namespace: GPArotation
Code
plot(rez_efa$n_fac)

Code
rez_efa$p

Final

We discarded:

  • Taste: Lone item + unstable
  • Affective_touch: Cross-loadings + unstable
  • Vomit: Less strongly associated
  • Itch: Less strongly associated, did not form consistent cluster
  • More controversial: Temp and Sex_arousal: similar cluster yet less reliable
data_allf <- select(data_all, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data1af <- select(data1a, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data1bf <- select(data1b, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data2f <- select(data2, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data3f <- select(data3, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data4f <- select(data4, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data5f <- select(data5, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data6f <- select(data6, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data7af <- select(data7a, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data7bf <- select(data7b, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data7cf <- select(data7c, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data8af <- select(data8a, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data8bf <- select(data8b, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data9f <- select(data9, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data10af <- select(data10a, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data10bf <- select(data10b, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data10cf <- select(data10c, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data11af <- select(data11a, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data11bf <- select(data11b, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
Code
rez_hclust <- make_hclust(data_allf)

plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
rez_hclust$p

Code
rez_ega <- make_ega(data_allf, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))

stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.528 0.042 0.017 0.041 0.026 0.042 0.134
Breathing 0.528 0.203 0.065 0.054 0.01  0.17  0.01 
Hungry 0.144 0.525 0.083 0.008 0.013 0.096 0.053
Thirsty 0.098 0.525 0.128 0.018 0.018 0.104 0.017
Urinate 0.032 0.192 0.572 0.043 0.025 0.116 0    
Defecate 0.069 0.072 0.572 0.087 0.062 0.115 0    
Sneeze 0.057 0.009 0.089 0.563 0.104 0.087 0.03 
Cough 0.054 0.021 0.036 0.563 0.176 0.12  0.078
Wind 0.014 0.024 0.071 0.117 0.576 0.074 0.067
Burp 0.032 0.015 0.016 0.179 0.576 0.098 0.127
Muscles 0.125 0.035 0.067 0.082 0.107 0.494 0.16 
Pain 0.054 0.135 0.09  0.066 0.008 0.494 0.198
Bruise 0.042 0     0     0.051 0.092 0.281 0.488
Blood_Sugar 0.075 0.057 0     0.023 0.035 0.065 0.488
Code
rez_ega$p

Code
rez_efa <- make_efa(data_allf, n = 3, cols = 
                      c("F1" = colors[["Wind"]], "F2" = colors[["Defecate"]], 
                        "F3" = colors[["Heart"]]))
plot(rez_efa$n_fac)

Code
rez_efa$p

Code
make_cfa <- function(data) {
  
  # Main models
  m_g <- '
IAS =~ Burp + Wind + Cough + Sneeze + Hungry + Thirsty  + Urinate + Defecate + Heart + Breathing + Pain + Muscles 
  '
  m_ega <- '
    CoughSneeze =~ Cough + Sneeze
    UrinateDefecate =~ Urinate + Defecate
    BruiseBloodsugar =~ Bruise + Blood_Sugar
    WindBurp =~ Wind + Burp
    HeartBreath =~ Heart + Breathing
    MusclesPain =~ Muscles + Pain
    HungryThirsty =~ Hungry + Thirsty
  '
  
  m_efa <- '
    UrinateDefecate =~ Urinate + Defecate
    Sudden =~ Cough + Sneeze + Wind + Burp
    Homeostasis =~ Heart + Breathing + Muscles + Pain + Hungry + Thirsty + Bruise + Blood_Sugar
  '
    
    
  m_hclust <- '
    HeartBreath =~ Heart + Breathing
    BruiseBloodsugar =~ Bruise + Blood_Sugar
    HungryThirsty =~ Hungry + Thirsty
    Sudden =~ Cough + Sneeze + Wind + Burp
    Unpleasant =~ Urinate + Defecate + Muscles + Pain
  '
  
  # Test bifactor models
  m_egag <- paste0(m_ega, "\nIAS =~ CoughSneeze + WindBurp + BruiseBloodsugar + UrinateDefecate + HeartBreath + MusclesPain")
  m_efag <- paste0(m_efa, "\nIAS =~ UrinateDefecate + Sudden + Homeostasis")
  m_hclustg <- paste0(m_hclust, "\nIAS =~ HeartBreath + BruiseBloodsugar + HungryThirsty + Sudden + Unpleasant")

  fit_g <- lavaan::cfa(m_g, data = data, std.lv = TRUE)
  fit_ega <- lavaan::cfa(m_ega, data = data, std.lv = TRUE)
  fit_egag <- lavaan::cfa(m_egag, data = data, std.lv = TRUE)
  fit_efa <- lavaan::cfa(m_efa, data = data, std.lv = TRUE)
  fit_efag <- lavaan::cfa(m_efag, data = data, std.lv = TRUE)
  fit_hclust <- lavaan::cfa(m_hclust, data = data, std.lv = TRUE)
  fit_hclustg <- lavaan::cfa(m_hclustg, data = data, std.lv = TRUE)
  
  # Test
  bf <- test_bf(fit_ega, fit_g, fit_egag, fit_efa, fit_efag, fit_hclust, fit_hclustg) |> 
    as.data.frame() |> 
    select(Name = Model, BF)
  
  # test <- lavaan::lavTestLRT(fit_g, fit_ega, fit_egag, fit_efa, fit_efag, fit_hclust, fit_hclustg)
  test <- compare_performance(fit_g, fit_ega, fit_egag, fit_efa, fit_efag, fit_hclust, fit_hclustg,
                      metrics = c("BIC", "RMSEA", "Chi2", "CFI")) |> 
    select(-BIC_wt, -Model) |> 
    merge(bf, by = "Name") |>
    arrange(BIC) |> 
    gt::gt() |> 
    gt::tab_header(title = "CFA Model Comparison") |>
    gt::data_color(
      columns = c(BIC, Chi2, RMSEA),
      method = "numeric",
      palette = c("green", "white")) |>
    gt::data_color(
      columns = c(CFI),
      method = "numeric",
      palette = c("white", "green")) |>
    gt::data_color(
      columns = c(BF),
      method = "numeric",
      palette = c("red", "yellow", "white", "yellow", "green"),
      domain = c(0, 0.3, 1, 3, 10000000)) |>
    gt::fmt_number(columns = c(RMSEA, Chi2, CFI), decimals = 3) 

  
  param  <- parameters::parameters(fit_ega)
  perf <- performance::performance(fit_ega)
  
  # Graph
  nodes <- tidySEM::get_nodes(fit_ega) |> 
    mutate(name = str_replace(name, "_", " "))
  edges <- tidySEM::get_edges(fit_ega) |> 
    mutate(from = str_replace(from, "_", " "),
           to = str_replace(to, "_", " "), 
           est_std = as.numeric(est_std)) |> 
    filter(from != to)
  
  g <- tidygraph::tbl_graph(nodes = nodes, edges = edges, directed = TRUE) 
  set.seed(1)
  layout <- ggraph::create_layout(g, layout = "fr", weights = est_std)
  
  p_cfa <- ggraph(layout) +
    geom_edge_bend(aes(edge_width = abs(est_std), edge_alpha = abs(est_std),
                       label = format_value(est_std, zap_small = TRUE, lead_zero = FALSE),
                       filter = arrow == "none"), 
                   color = "forestgreen",
                   angle_calc = "along", 
                   check_overlap = TRUE,
                   label_dodge = unit(0, "mm"),
                   label_size = unit(3, "mm"),
                   strength = 0.3) +
    geom_edge_link(aes(edge_width = abs(est_std), edge_color = abs(est_std), 
                       label = format_value(est_std, zap_small = TRUE, lead_zero = FALSE),
                      filter = arrow != "none"), 
                   angle_calc = "along", 
                   label_dodge = unit(3, "mm"),
                   arrow = arrow(length = unit(2.5, "mm"), type = "open", ends = "first"), 
                   lineend = "round",
                   linejoin = "bevel",
                   start_cap = circle(9, 'mm')) +
    geom_node_point(aes(color = name, shape=shape), alpha = 0.97, size = 20) +
    geom_node_text(data = filter(layout, shape == "rect"), aes(label = name), 
                   size = 3, color = "white", fontface = "bold") +
    scale_size_continuous(range = c(10, 20)) +
    scale_edge_width(range = c(0.3, 4)) +
    scale_edge_alpha(range = c(0.1, 0.9)) +
    scale_edge_color_gradient(low = "#EEEEEE", high = "#212121") +
    scale_color_manual(values = c(colors, 
                               "CoughSneeze" = "#CDDC39", 
                               "UrinateDefecate" = "#795548", "BruiseBloodsugar" = "#673AB7",
                               "WindBurp" = "#009688", "HeartBreath" = "#F44336",
                               "MusclesPain" = "#2196F3", "HungryThirsty" = "#FF9800")) +
    scale_shape_manual(values = c("oval" = 16, "rect" = 15)) +
    labs(title = "Confirmatory Factor Analysis (CFA)", subtitle = "Method = Maximum Likelihood") +
    theme_void() +
    theme(legend.position = "none",
          plot.title = element_text(face="bold"))
  p_cfa
  
  list(test = test, param = param, perf = perf, p=p_cfa)
}


rez_cfa0 <- make_cfa(data_allf)
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: When comparing models, please note that probably not all models were fit
  from same data.
When comparing models, please note that probably not all models were fit
  from same data.
Code
rez_cfa0$test
CFA Model Comparison
Name BIC RMSEA Chi2 CFI BF
fit_ega -158722.7 0.036 2,290.274 0.984 NA
fit_egag -154419.8 0.055 6,738.245 0.952 0
fit_hclust -148000.0 0.079 13,126.977 0.905 0
fit_hclustg -146829.7 0.079 14,349.064 0.896 0
fit_g -145655.5 0.125 26,541.561 0.783 0
fit_efa -144933.2 0.083 16,266.333 0.883 0
fit_efag -144933.2 0.083 16,266.333 0.883 0
Code
rez_cfa0$perf
# Indices of model performance

Chi2(56) | p (Chi2) | Baseline(91) | p (Baseline) |   GFI |  AGFI |   NFI
-------------------------------------------------------------------------
2290.274 |   < .001 |    1.380e+05 |       < .001 | 0.990 | 0.981 | 0.983

Chi2(56) |  NNFI |   CFI | RMSEA |      RMSEA  CI | p (RMSEA) |   RMR |  SRMR
-----------------------------------------------------------------------------
2290.274 | 0.974 | 0.984 | 0.036 | [0.034, 0.037] |    > .999 | 0.001 | 0.020

Chi2(56) |   RFI |  PNFI |   IFI |   RNI | Loglikelihood |        AIC
---------------------------------------------------------------------
2290.274 | 0.973 | 0.605 | 0.984 | 0.984 |     79615.220 | -1.591e+05

Chi2(56) |        BIC | BIC_adjusted
------------------------------------
2290.274 | -1.587e+05 |   -1.589e+05
Code
rez_cfa0$param
# Loading

Link                            | Coefficient |       SE |       95% CI
-----------------------------------------------------------------------
CoughSneeze =~ Cough            |        0.17 | 1.16e-03 | [0.17, 0.17]
CoughSneeze =~ Sneeze           |        0.15 | 1.10e-03 | [0.15, 0.15]
UrinateDefecate =~ Urinate      |        0.17 | 1.23e-03 | [0.17, 0.17]
UrinateDefecate =~ Defecate     |        0.16 | 1.19e-03 | [0.16, 0.16]
BruiseBloodsugar =~ Bruise      |        0.22 | 2.11e-03 | [0.21, 0.22]
BruiseBloodsugar =~ Blood_Sugar |        0.16 | 1.94e-03 | [0.15, 0.16]
WindBurp =~ Wind                |        0.18 | 1.29e-03 | [0.18, 0.18]
WindBurp =~ Burp                |        0.20 | 1.32e-03 | [0.20, 0.20]
HeartBreath =~ Heart            |        0.14 | 1.49e-03 | [0.14, 0.14]
HeartBreath =~ Breathing        |        0.17 | 1.44e-03 | [0.17, 0.18]
MusclesPain =~ Muscles          |        0.15 | 1.19e-03 | [0.15, 0.15]
MusclesPain =~ Pain             |        0.15 | 1.23e-03 | [0.14, 0.15]
HungryThirsty =~ Hungry         |        0.18 | 1.62e-03 | [0.18, 0.18]
HungryThirsty =~ Thirsty        |        0.18 | 1.58e-03 | [0.18, 0.19]

Link                            |      z |      p
-------------------------------------------------
CoughSneeze =~ Cough            | 147.35 | < .001
CoughSneeze =~ Sneeze           | 136.55 | < .001
UrinateDefecate =~ Urinate      | 135.80 | < .001
UrinateDefecate =~ Defecate     | 136.25 | < .001
BruiseBloodsugar =~ Bruise      | 102.29 | < .001
BruiseBloodsugar =~ Blood_Sugar |  80.90 | < .001
WindBurp =~ Wind                | 141.15 | < .001
WindBurp =~ Burp                | 150.94 | < .001
HeartBreath =~ Heart            |  95.60 | < .001
HeartBreath =~ Breathing        | 119.51 | < .001
MusclesPain =~ Muscles          | 124.42 | < .001
MusclesPain =~ Pain             | 119.27 | < .001
HungryThirsty =~ Hungry         | 112.05 | < .001
HungryThirsty =~ Thirsty        | 115.23 | < .001

# Correlation

Link                                | Coefficient |       SE |       95% CI
---------------------------------------------------------------------------
CoughSneeze ~~ UrinateDefecate      |        0.56 | 5.94e-03 | [0.55, 0.57]
CoughSneeze ~~ BruiseBloodsugar     |        0.55 | 7.32e-03 | [0.53, 0.56]
CoughSneeze ~~ WindBurp             |        0.72 | 4.80e-03 | [0.71, 0.73]
CoughSneeze ~~ HeartBreath          |        0.52 | 6.73e-03 | [0.51, 0.53]
CoughSneeze ~~ MusclesPain          |        0.67 | 5.99e-03 | [0.65, 0.68]
CoughSneeze ~~ HungryThirsty        |        0.46 | 7.02e-03 | [0.45, 0.47]
UrinateDefecate ~~ BruiseBloodsugar |        0.40 | 7.83e-03 | [0.38, 0.41]
UrinateDefecate ~~ WindBurp         |        0.51 | 6.10e-03 | [0.50, 0.52]
UrinateDefecate ~~ HeartBreath      |        0.54 | 6.71e-03 | [0.52, 0.55]
UrinateDefecate ~~ MusclesPain      |        0.65 | 6.14e-03 | [0.64, 0.66]
UrinateDefecate ~~ HungryThirsty    |        0.65 | 6.15e-03 | [0.64, 0.66]
BruiseBloodsugar ~~ WindBurp        |        0.58 | 7.09e-03 | [0.57, 0.59]
BruiseBloodsugar ~~ HeartBreath     |        0.50 | 8.08e-03 | [0.48, 0.51]
BruiseBloodsugar ~~ MusclesPain     |        0.74 | 7.37e-03 | [0.73, 0.76]
BruiseBloodsugar ~~ HungryThirsty   |        0.45 | 8.30e-03 | [0.43, 0.47]
WindBurp ~~ HeartBreath             |        0.46 | 6.86e-03 | [0.44, 0.47]
WindBurp ~~ MusclesPain             |        0.63 | 6.03e-03 | [0.62, 0.65]
WindBurp ~~ HungryThirsty           |        0.43 | 7.01e-03 | [0.42, 0.45]
HeartBreath ~~ MusclesPain          |        0.64 | 6.88e-03 | [0.63, 0.66]
HeartBreath ~~ HungryThirsty        |        0.64 | 6.90e-03 | [0.63, 0.66]
MusclesPain ~~ HungryThirsty        |        0.63 | 6.90e-03 | [0.62, 0.65]

Link                                |      z |      p
-----------------------------------------------------
CoughSneeze ~~ UrinateDefecate      |  94.23 | < .001
CoughSneeze ~~ BruiseBloodsugar     |  74.98 | < .001
CoughSneeze ~~ WindBurp             | 150.87 | < .001
CoughSneeze ~~ HeartBreath          |  77.12 | < .001
CoughSneeze ~~ MusclesPain          | 111.01 | < .001
CoughSneeze ~~ HungryThirsty        |  65.50 | < .001
UrinateDefecate ~~ BruiseBloodsugar |  51.02 | < .001
UrinateDefecate ~~ WindBurp         |  83.50 | < .001
UrinateDefecate ~~ HeartBreath      |  79.79 | < .001
UrinateDefecate ~~ MusclesPain      | 105.72 | < .001
UrinateDefecate ~~ HungryThirsty    | 105.96 | < .001
BruiseBloodsugar ~~ WindBurp        |  81.97 | < .001
BruiseBloodsugar ~~ HeartBreath     |  61.56 | < .001
BruiseBloodsugar ~~ MusclesPain     | 100.55 | < .001
BruiseBloodsugar ~~ HungryThirsty   |  54.31 | < .001
WindBurp ~~ HeartBreath             |  66.62 | < .001
WindBurp ~~ MusclesPain             | 105.09 | < .001
WindBurp ~~ HungryThirsty           |  61.95 | < .001
HeartBreath ~~ MusclesPain          |  93.45 | < .001
HeartBreath ~~ HungryThirsty        |  93.47 | < .001
MusclesPain ~~ HungryThirsty        |  91.62 | < .001
Code
rez_cfa0$p

Figure

The structural analysis seem to converge on the idea of small clusters (of pairs or triplets) that are potentially inter-related parts of larger clusters (although with unstable associations). Best way to assess the values of this new granular structure is to assess its predictive value against convergent measures, and see if it’s superior to a unique score (-> Study 2).

Code
fig2 <- (rez_hclust$p | rez_ega$p) / (rez_efa$p | rez_cfa0$p) 
ggsave("figures/Figure2.png", fig2, width=14, height=14, dpi=300, bg="white")

Structure Invariance

Initial

Code
rez_hclust <- make_hclust(data1a)
rez_ega <- make_ega(data1a, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data1a, n = 4, cols = c("F1" = colors[["Wind"]], "F2" = colors[["Defecate"]], 
                                            "F3" = colors[["Heart"]], "F4" = colors[["Itch"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7 8
Heart 0.504 0.125  0.046 0     0.093 0.114 0      0    
Breathing 0.504 0.14   0.032 0.098 0.158 0.123 0      0.036
Hungry 0.2   0.511  0.12  0.004 0.028 0.064 0      0.034
Thirsty 0.074 0.511  0.166 0     0.068 0.13  0.005  0.053
Defecate 0.044 0.093  0.632 0.107 0.077 0.036 0.059 −0.045
Urinate 0.026 0.253  0.421 0.035 0.167 0.035 0.046 −0.02 
Taste 0.032 0.018  0.21  0.105 0.18  0.149 0.082  0.028
Sneeze 0     0.006  0.078 0.676 0.094 0.001 0.092  0.044
Cough 0.053 0      0.014 0.377 0.11  0     0.171  0.121
Vomit 0.088 0      0.177 0.299 0.066 0.067 0      0.054
Muscles 0.084 0.022  0.13  0.03  0.454 0.032 0.029  0.093
Pain 0.126 0.032  0.133 0.067 0.438 0.195 0      0.143
Temp 0.086 0.056  0.117 0.125 0.264 0.245 0      0.038
Sex_arousal 0.103 0.042  0.055 0.006 0.069 0.435 0.119  0.007
Affective_touch 0.059 0.086  0.058 0.027 0.205 0.435 0.011  0.152
Wind 0     0.007  0.189 0.051 0     0.261 0.571  0.058
Burp 0     0      0.007 0.203 0.034 0     0.571  0.105
Blood_Sugar 0.013 0     −0.02  0.04  0.033 0.033 0.015  0.455
Bruise 0     0.108 −0.044 0.08  0.193 0.044 0.068  0.401
Itch 0.034 0      0.027 0.076 0.074 0.224 0.07   0.388
Code
plot(rez_efa$n_fac)

Final

Code
rez_hclust <- make_hclust(data1af)
rez_ega <- make_ega(data1af, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data1af, n = 4, cols = c("F1" = colors[["Wind"]], "F2" = colors[["Defecate"]], 
                                            "F3" = colors[["Heart"]], "F4" = colors[["Itch"]]))
rez_cfa1a  <- make_cfa(data1af)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.511 0.128  0.036 0.009 0.01  0.108  0    
Breathing 0.511 0.16   0.048 0.068 0.015 0.151  0.045
Hungry 0.202 0.514  0.105 0.015 0     0.042  0.048
Thirsty 0.091 0.514  0.178 0     0.018 0.041  0.089
Urinate 0.029 0.253  0.565 0.038 0.05  0.185 −0.017
Defecate 0.08  0.108  0.565 0.098 0.087 0.09  −0.042
Sneeze 0.033 0.02   0.12  0.574 0.117 0.091  0    
Cough 0.071 0      0.021 0.574 0.171 0.06   0.157
Wind 0.021 0.024  0.136 0.071 0.577 0.028  0.089
Burp 0.012 0      0.008 0.22  0.577 0.041  0.089
Muscles 0.091 0.042  0.118 0.058 0.04  0.523  0.125
Pain 0.183 0.044  0.107 0.061 0.014 0.523  0.18 
Bruise 0.01  0.128 −0.031 0.059 0.093 0.222  0.502
Blood_Sugar 0.033 0     −0.013 0.052 0.032 0.053  0.502
Code
plot(rez_efa$n_fac)

Code
rez_cfa1a$test
CFA Model Comparison
Name BIC RMSEA Chi2 CFI BF
fit_ega -1635.181 0.026 72.825 0.988 NA
fit_egag -1613.162 0.059 180.403 0.924 1.654961e-05
fit_hclustg -1594.961 0.065 210.828 0.905 1.846821e-09
fit_hclust -1592.739 0.062 182.492 0.921 6.080789e-10
fit_efag -1564.540 0.073 253.472 0.877 4.576430e-16
fit_efa -1564.540 0.073 253.472 0.877 4.576430e-16
fit_g -1538.631 0.106 325.704 0.784 1.082552e-21

TODO.

Model Performance

Code
rbind(
  mutate(rez_cfa0$perf, Sample = "All Samples"),
  mutate(rez_cfa1a$perf, Sample = "Sample 1a")
  # TODO
) |> 
  select(Sample, Chi2, RMSEA, CFI, "SRMR") 
# Indices of model performance

Sample      |     Chi2 | RMSEA |   CFI |  SRMR
----------------------------------------------
All Samples | 2290.274 | 0.036 | 0.984 | 0.020
Sample 1a   |   72.825 | 0.026 | 0.988 | 0.029

Scores

Note the usage of descriptive factor names relating directly to the items to avoid abstraction.

Code
# Add empirical variables
add_empirical <- function(data, sample = "Sample 1a") {
  x <- data.frame(
    Original = rowMeans(data),
    HungryThirsty = (data$Hungry + data$Thirsty) / 2,
    MusclesPain = (data$Muscles + data$Pain) / 2,
    WindBurp = (data$Wind + data$Burp) / 2,
    UrinateDefecate = (data$Urinate + data$Defecate) / 2,
    BreathingHeart = (data$Heart + data$Breathing) / 2)
  
  if("Blood_Sugar" %in% names(data)) {
    x$BruiseBlood <- (data$Blood_Sugar + data$Bruise) / 2
  } else {
    x$BruiseBlood <- data$Bruise
  }
  
  if("Cough" %in% names(data)) {
    x$CoughSneeze <- (data$Cough + data$Sneeze) / 2
  } else {
    x$CoughSneeze <- data$Sneeze
  }
  x$Sample <- sample
  x
}

scores <- list(
  sample1a = add_empirical(data1a, "Sample 1a"),
  sample1b = add_empirical(data1b, "Sample 1b"),
  sample2 = add_empirical(data2, "Sample 2"),
  sample3 = add_empirical(data3, "Sample 3"),
  sample4 = add_empirical(data4, "Sample 4"),
  sample5 = add_empirical(data5, "Sample 5"),
  sample6 = add_empirical(data6, "Sample 6"),
  sample7a = add_empirical(data7a, "Sample 7a"),
  sample7b = add_empirical(data7b, "Sample 7b"),
  sample7c = add_empirical(data7c, "Sample 7c"),
  sample8a = add_empirical(data8a, "Sample 8a"),
  sample8b = add_empirical(data8b, "Sample 8b"),
  sample9 = add_empirical(data9, "Sample 9"),
  sample10a = add_empirical(data10a, "Sample 10a"),
  sample10b = add_empirical(data10b, "Sample 10b"),
  sample10c = add_empirical(data10c, "Sample 10c"),
  sample11a = add_empirical(data11a, "Sample 11a"),
  sample11b = add_empirical(data11b, "Sample 11b")
)

save(scores, file = "../data/scores.RData")